        # rm(list=ls()) ; gc()


      #Article:    Legislature Resizing with Rent-Seeking Politicians: The Impact of Executive-Legislature Coalitions			
      #Author:     Anderson Frey <anderson.frey@rochester.edu>			    

      # This script produces all Tables and Figures in the Article
      # Tables are produced in a way that they can be directly copied/pasted in a LaTex document
      
      # It needs the following files:  data_main.csv  | data_reelection.csv 

      # Please upload all libraries and functions below before attempting to produce any Tables or Figures

      #Set the working directory to the folder that contains the replication files, use the function below
      
      #setwd()
        
      # Libraries ========================================================================================================
      
        
        library(data.table) 
        library(dplyr) 
        library(lfe)
        library(rdrobust)
        library(ggplot2) 
        library(rddensity)
        library(gridExtra)
        library(fastDummies)
        library(tidyr)
        library(boot)
        library(ggpubr)
        
        
        
      #===================================================================================================================
      
        
      # Functions ========================================================================================================
      
  
        # Regression output in LaTex format
        .LatexTable <- function(data,rownames,spaces){
          cols <- ncol(data)*2+1
          res <- data.frame(matrix(0,nrow(data),cols))
          res[,] <- "&"
          fill <- seq(from=2,to=cols,by=2)
          res[,fill] <- data
          res <- cbind(rownames,res)
          lastone <- rep("\tabularnewline",nrow(res))
          lastone[spaces] <- "\tabularnewline[0.2cm]"
          lastone[nrow(res)] <- "\tabularnewline[0.2cm]"
          res[,ncol(res)] <- lastone
          names(res) <- NULL
          return(res)
        }
        
        # Plot of marginal effects
        .NewPlot <- function(data,r,v,xvar,minD,maxD,bins1,bins2,shift,bb,x.title,y.title,
                             p.title,degree,multi,oneside,bb3, pt){
          
          
          Di <- matrix(0,(bins1),3)
          Di[,1] <- seq(from=minD, to=maxD, length.out=(bins1))
          label <- paste("t",xvar,sep=":")
          
          for (i in 1:nrow(Di)){
            Di[i,2] <- r["t",1] + r[label,1]*Di[i,1] 
            Di[i,3] <-(  v["t","t"] + (Di[i,1]**2)*v[label,label] + 2*Di[i,1]*v["t",label] )^.5
          }
          
          MT <- data.frame(Di) ; names(MT) <- c("D","coeff","sd")
          MT$h <- MT$coeff+1.96*MT$sd
          MT$l <- MT$coeff-1.96*MT$sd  
          
          jazz <- data.frame(matrix(0,(bins2),2))
          xx <- seq(from=min(MT$D), to=(max(MT$D)), length.out=(bins2+1))
          
          data <- data.frame(data)
          data$x <- data[,match(xvar,names(data))]
          
          i <- 2
          
          for (i in 1:nrow(jazz)){
            top <- xx[i+1]
            bot <- xx[i]
            if(is.na(top))  top <- xx[i+1] + 1
            jj <- subset(data, x >= bot & x < top )
            if(oneside==1) jj <- subset(data, x >= bot & x < top & t == 1)
            if(oneside==1) {datax <- subset(data,  t == 1)}else{datax <- data}
            jazz[i,1] <- (xx[i]+xx[i+1])/2
            
            jazz[i,2] <- nrow(jj)/nrow(datax)
          }
          names(jazz) <- c("D","dens")
          
          bb2 <- bb+shift
          
          if(pt==1){
            ggplot(data=MT, aes(x=D, y=coeff, ymin=0))   +    theme_bw() +
              theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
              geom_line(data=MT, aes(x=D, y=coeff+shift),colour="black",  linewidth=1, alpha=.75) +
              theme(plot.title = element_text(family = "A", size=13))+
              geom_hline(yintercept=(0+shift), lty=3)+
              xlab(x.title) + ylab(y.title) +
              theme(axis.text = element_text(colour="black", size=13, family="A")) +
              theme(axis.title = element_text(colour="black", size=13, family="A")) +
              geom_line(data=MT, aes(x=D, y=h+shift),colour="black",  linetype="dashed",linewidth=.75, alpha=.75) +
              geom_line(data=MT, aes(x=D, y=l+shift),colour="black",  linetype="dashed",linewidth=0.75, alpha=.75) +
              ggtitle(p.title) +
              geom_bar(data=jazz, aes(x=D, y=dens*multi, group = 1), stat="identity", position="dodge",
                       size =0.1,color="gray35",fill="gray85", alpha=0.3) +
              scale_y_continuous(breaks=bb2, labels=bb )+
              scale_x_continuous(breaks=bb3)+
              coord_cartesian(ylim=c(bb2[1],bb2[length(bb2)]))
          }else{
            ggplot(data=MT, aes(x=D, y=coeff, ymin=0))   +    theme_bw() +
              theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
              geom_line(data=MT, aes(x=D, y=coeff+shift),colour="black",  linewidth=1, alpha=.75) +
              theme(plot.title = element_text(family = "A", size=13))+
              geom_hline(yintercept=(0+shift), lty=3)+
              xlab(x.title) + ylab(y.title) +
              theme(axis.text = element_text(colour="black", size=13, family="A")) +
              theme(axis.title = element_text(colour="black", size=13, family="A")) +
              geom_line(data=MT, aes(x=D, y=h+shift),colour="black",  linetype="dashed",linewidth=.75, alpha=.75) +
              geom_line(data=MT, aes(x=D, y=l+shift),colour="black",  linetype="dashed",linewidth=0.75, alpha=.75) +
              geom_bar(data=jazz, aes(x=D, y=dens*multi, group = 1), stat="identity", position="dodge",
                       size =0.1,color="gray35",fill="gray85", alpha=0.3) +
              scale_y_continuous(breaks=bb2, labels=bb )+
              scale_x_continuous(breaks=bb3)+
              coord_cartesian(ylim=c(bb2[1],bb2[length(bb2)]))
          }
          
        }
        
        # For bootstrapping difference in coefficients from 2 RDD regressions
        .DiffRDD <- function( base, indices, b1,b2, controls, variable){
          
          library(fastDummies)
          library(rdrobust)
          library(dplyr)
          library(data.table)
          
          bb <- data.table( base[indices,] )
          ddd <- bb %>% select((fe))
          datafex <- dummy_cols(ddd)
          datafex <- datafex[,3:ncol(datafex)]
          
          
          bb$i <- bb %>% select(all_of(variable))
          
          bb$ti <- bb$t*bb$i
          bb$tir <- bb$t*bb$i*bb$run
          bb$ir <- bb$i*bb$run
          dataconx <- bb %>% select(all_of(controls))
          cv <- cbind(bb$i,bb$ti,bb$tir,bb$ir,dataconx,datafex)
          
          reg <- rdrobust(bb$var,bb$run,  b=b2, h=b1,  covs=cv,  masspoints="off") 
          a <- reg$coef[1]
          bb$i <- 1-bb$i
          
          bb$ti <- bb$t*bb$i
          bb$tir <- bb$t*bb$i*bb$run
          bb$ir <- bb$i*bb$run
          
          cv <- cbind(bb$i,bb$ti,bb$tir,bb$ir,dataconx,datafex)
          
          reg <- rdrobust(bb$var,bb$run,  b=b2, h=b1, covs=cv, masspoints="off") 
          b <- reg$coef[1]
          
          return(b-a)
          
        }
        
        
      #===================================================================================================================
      
      
      # Setup ============================================================================================================
        
        # set fonts for plots
        windowsFonts( A=windowsFont("Liberation Sans")) # Font for ggplot2 
        
        
        # Open file
        base <- read.csv("data_main.csv")
   
        
        # Adjust covariates
        base$lpoverty <- log(base$poverty)
        base$lvoters <- log(base$present08)
        base$lArea <- log(base$Area)
        base$rent.inc <- 0
        base$rent.inc[base$m.party.last %in% c(15,11,22,14)] <- 1
        base$pt_base <- ifelse(base$m.party %in% c(11,12,13,14,15,22,40),1,0)
        base$v <- (base$vote_win*base$last.w + base$op.vote_win*base$last.l )*100 / base$tv_winner 
        base$run <- base$residual*100 / base$tv_winner
        base$run <- (2*base$last.w-1)*base$run
        
        # Define treatment dummy and fixed-effects
        base$t <- base$last.w
        base$fe <- base$seats08
        
        # Parameters for optimal bandwidth algorithm 
        btype <- "cercomb1"
        qq <- 5
        vart <- "hc3"
        sdd <- sd(abs(base$run))

        
        # List of covariates
        controls <- c("qe","lgdppc","lbpc","lvoters","lpoverty","rent.inc","lArea","pt_base")
    
       
        
      #===================================================================================================================

    
      # Figure 1 ==========================================================================================================
          
          
          m <- fread("data_reelection.csv")
          m$year <- ifelse(m$election==2012,1,0)
          
          xx1 <- group_by(m,year,change) %>% summarise(run=mean(run)) 
          xx1$share <- xx1$run*100
          xx2 <- group_by(m,year,change, run) %>% summarise(win=mean(win)) 
          xx2 <- filter(xx2, run==1)
          xx2$share <- xx2$win*100
          xx2$run <- NULL
          
          xx1$type <- "Run Again"
          xx2$type <- "Re-elected"
          
          plot_data <- rbind(xx1,xx2)
          plot_data$x <- "Attempted Reelection \n in 2008 (Placebo)"
          plot_data$x[plot_data$year == 1] <- "Attempted Reelection \n in 2012"
          plot_data$x[plot_data$year == 0 & plot_data$type == "Re-elected"] <- "Reelected \n in 2008 (Placebo)"
          plot_data$x[plot_data$year == 1 & plot_data$type == "Re-elected"] <- "Reelected \n in 2012"
          plot_data$g <- as.factor(plot_data$change)
          
          reg1 <- felm(run ~ (change) |  0 | 0 | reg, data=m[m$year == 0] ) ; summary(reg1)
          reg2 <- felm(run ~ (change) |  0 | 0 | reg, data=m[m$year == 1] ) ; summary(reg2)
          reg3 <- felm(win ~ (change) |  0 | 0 | reg, data=m[m$year == 0 & m$run==1] ) ; summary(reg3)
          reg4 <- felm(win ~ (change) |  0 | 0 | reg, data=m[m$year == 1 & m$run==1] ) ; summary(reg4)
          
          label1 <- paste0("Difference: ",round(summary(reg1)$coefficients["change",1]*100,1),"%")
          label2 <- paste0("Difference: ",format(round(summary(reg2)$coefficients["change",1]*100,1),nsmall=1),"%*")
          label3 <- paste0("Difference: ",format(round(summary(reg3)$coefficients["change",1]*100,1),nsmall=1),"%")
          label4 <- paste0("Difference: ",format(round(summary(reg4)$coefficients["change",1]*100,1),nsmall=1),"%*")
          
          sz <- 12
          alphayr <- c("0.5","0.5","1","1","0.5","0.5","1","1")
          
          ggplot(data=plot_data, aes( x=x, y=share, group=g))   +   theme_bw() +
            ylim(0,100)+
            theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
            geom_bar(data=plot_data, aes( x=x,y=share, group=g, fill=g, alpha=factor(alphayr)), 
                     stat="identity", position_dodge()) +
            xlab("Share of Parties in a Coincident Coalition") + ylab("Percentage of Councilors") +
            theme(axis.text.y = element_text(colour="black", size=sz, family="A")) +
            theme(axis.title.x = element_blank()) +
            theme(axis.title = element_text(colour="black", size=sz, family="A")) +
            theme(legend.title = element_text(colour="black", size=sz, family="A")) +
            scale_fill_manual(values=c("gray25","gray75"),name="Group of Municipalities",
                              labels=c("Council size remained unchanged","Council size was increased"))+
            theme(legend.position = c(0.8,0.85))+
            scale_alpha_manual(values = c("0.5"=0.25, "1"=1), guide='none')+
            annotate(geom="text", x=1, y=85, label=label1,color="gray50",size=4, family="A")+
            annotate(geom="text", x=2, y=85, label=label2,color="black",size=4, family="A") +
            annotate(geom="text", x=3, y=65, label=label3,color="gray50",size=4, family="A")+
            annotate(geom="text", x=4, y=65, label=label4,color="black",size=4, family="A")
          
          rm(m,xx1,xx2,reg1,reg2,reg3,reg4,label1,label2,label3,label4,plot_data)
        
          
      #===================================================================================================================

        
      # Table 1 ==========================================================================================================
        
          
        # keeps only observations where the close election was: coalition vs. opposition
        new <- filter(base, MainSample==1) 
       
        # runs the regressions
        for (out in 1:2){
          
          if(out==1){new$var <- new$outcome}else{new$var <- new$t.mayor}
          res <- data.frame(matrix(0,5,3))
          
          ddd <- new %>% select(fe)
          datafe <- dummy_cols(ddd)
          datafe <- datafe[,3:ncol(datafe)]
          datacon <- new %>% select(all_of(controls))
          ddd <- new %>% select(reg)
          datareg <- dummy_cols(ddd)
          datareg <- datareg[,3:ncol(datareg)]

          b <- rdbwselect(new$var, new$run, p=1, q=qq, bwselect = btype)$bws ; b/sdd
          b1 <- b[1]
          b2 <- b1*(b[3]/b[1])
          
          
          for (i in 1:3){
            
            
            cv <- datafe
            if(i == 2) cv <- cbind(datafe,datacon)
            if(i == 3) cv <- cbind(datafe,datacon,datareg)
            
            reg0 <- rdrobust(new$var,new$run, b=b2, h=b1, vce=vart , q=qq) 
            reg <- rdrobust(new$var,new$run,  b=b2, h=b1, covs=cv, q=qq, vce=vart) #; summary(reg) 
            reg2 <- rdrobust(new$var,new$run, covs=cv, b=b2, h=b1, q=qq, level=c(90), vce=vart) 
            co <- round(reg$coef[1],3)
            se <- round(reg$se[3],3)
            ci1 <- round(reg$ci[3,1],2)
            ci2 <- round(reg$ci[3,2],2)
            int <- round(reg0$tau_cl[1],3)
            sam <- sum(reg$N_h)
            
            stars <- ""
            if(reg2$ci[3,1]*reg2$ci[3,2]>0) stars <- "+"
            if(reg$ci[3,1]*reg$ci[3,2]>0) stars <- "*"
            
            res[1,i] <- paste(format(co,nsmall=3),stars,sep="")
            res[2,i] <- paste( "(" ,format(se,nsmall=3), ")",sep="")
            res[3,i] <- paste( "[" ,format(ci1,nsmall=2),",",format(ci2,nsmall=2), "]",sep="")
            res[4,i] <- format(round(b1/sdd,2),nsmall=2)
            res[5,i] <- format(sam,nsmall=0)
          }
          
          if(out==1) res1 <- res
          
        }
        
        # produces the LaTex Table
        resy <- cbind(res1,res)
        spaces <- c(2,3,4,5)
        nam <- c("RD Effect"," ","Robust C.I.","Bandwidth","Observations")
        print(.LatexTable(resy,nam,spaces), row.names=F)
        
        
      #===================================================================================================================

            
      # Figure 2 =========================================================================================================
        
        # keeps only observations where the close election was: coalition vs. opposition
        basex <- filter(base, MainSample==1)    
        
        # Keep only data with donor information
        new <- filter(basex,  !is.na(share)  )  
        new$var <- new$outcome
        
        new$patr <- ifelse(new$share > median(new$share),1,0)
    
        # Optimal bandwidth
        b <- rdbwselect(new$var, new$run, p=1, q=qq, bwselect = btype)$bws ; b
        b1 <- b[1]
        b2 <- b1*(b[3]/b[1])
        
        res <- data.frame(matrix(0,3,7))
        bb <- new
        
        # Calculate individual regressions for each subsample based on patronage
        for (i in 0:1){
          
          ddd <- bb %>% select(fe)
          datafex <- dummy_cols(ddd)
          datafex <- datafex[,3:ncol(datafex)]
          dataconx <- new %>% select(all_of(controls))
          ddd <- new %>% select(reg)
          datastax <- dummy_cols(ddd)
          datastax <- datastax[,3:ncol(datastax)]

          bb$i <- abs(i-bb$patr)
          bb$tr <- bb$t*bb$run
          bb$ti <- bb$t*bb$i
          bb$tir <- bb$t*bb$i*bb$run
          bb$ir <- bb$i*bb$run
          cv <- cbind(bb$i,bb$ir,bb$tir,bb$ti,datafex,dataconx,datastax)
          
          reg <- rdrobust(bb$var,bb$run,  b=b2, h=b1,covs=cv ,q=qq , masspoints="off") 
          reg2 <- rdrobust(bb$var,bb$run,  b=b2, h=b1,covs=cv , q=qq , level=90 , masspoints="off") 
     
          stars <- ""
          if(reg2$ci[3,2]*reg2$ci[3,1]>0) stars <- "+"
          if(reg$ci[3,2]*reg$ci[3,1]>0) stars <- "*"
          
          
          res[1+i,1] <- reg$coef[1]
          res[1+i,2] <- reg$ci[3,1]
          res[1+i,3] <- reg$ci[3,2]
          res[1+i,4] <- reg2$ci[3,1]
          res[1+i,5] <- reg2$ci[3,2]
          res[1+i,6] <- stars
          res[1+i,7] <- b1
          
        }
        
        res[3,1] <- res[2,1] - res[1,1]
        
        # Bootstrap the CI for the difference in coefficients
        set.seed(14618) 
        results <- boot(data=new,statistic=.DiffRDD , R=500,
                        b1=b1,b2=b2, 
                        controls=controls,variable="patr",
                        parallel = "snow", ncpus = 6)
        
        ci = boot.ci(results, type=c("perc","basic"), conf=.95) ;  ci 
        ci2 = boot.ci(results, type=c("perc","basic"), conf=.90) ;  ci2
        res[3,2] <-    ci$basic[4]
        res[3,3] <-    ci$basic[5]
        stars <- ""
        if(ci2$basic[4]*ci2$basic[5]>0) stars <- "+"
        if(ci$basic[4]*ci$basic[5]>0) stars <- "*"
        res[3,6] <- stars
        names(res) <- c("co","lo","hi","lo2","hi2","stars","band")
        res$nm <- c("Low \n Patronage","High \n Patronage","Difference")
        res$x <- as.factor(1:3)
        res_patr <- res
        
        # Left-side plot
        ph1 <- ggplot(data=res_patr, aes(x=x, y=co))   + theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
          geom_hline(yintercept=0, lty=3)+
          geom_point( size=5)+
          geom_errorbar(aes(ymin=lo,ymax=hi),width=0.0, linewidth=1.25)+
          scale_x_discrete(labels=res_patr$nm) +
          theme(axis.text = element_text(colour="black", size=13, family="A")) +
          theme(axis.title = element_blank())+
          ylab("Coefficient (RD Effect)") 
     
     
        
        # Restrict the data to observations within the optimal band
        data_short <- filter(new, abs(run) < b1)
        data_short$var <- data_short$outcome
        data_short$wght <- b1 - abs(data_short$run)
        
        formula2 <- as.formula( paste0( "var ~ t*run*share+", paste(controls,collapse="+") ,"|reg+fe|0|0"  ))
        reg <- felm( formula2,  weights=data_short$wght, data=data_short) ; summary(reg, robust=TRUE)
        
        
        rr <- summary(reg)$coefficients
        vv <- vcov(reg)
        
        vec <- round(seq(from= -2, to = 1,length.out=4),2)
        vec2 <- round(seq(from= -1, to = 1,length.out=5),1) #c(9,11,13,15,17,19,21)
        
        # Right-side plot
        ph2 <- .NewPlot(data=data_short,r=rr,v=vv,xvar="share",minD=-1,maxD=1,bins1=12,bins2=14,
                        shift=2,bb=vec,degree=1,multi=6,oneside=0,
                        x.title="Patronage Level",y.title="",p.title="",bb3=vec2,pt=0) 
    
        
        
        grid.arrange(ph1,ph2,ncol=2,widths=c(0.4, 0.6),left = text_grob("RDD Effect", size = 13,  family="A", rot = 90)) 
        
        rm(bb,new,rr,vv,data_short,res_patr,res)
     
      #===================================================================================================================
        
        
      
      
      

      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
